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Abstract. We consider the radial transport of test particles due to the E x B 
drift motion in the guiding center approximation. Using an explicit expression 
to modify the electrostatic potential, we show that it is possible to construct a 
transport barrier which suppresses radial transport. We propose an algorithm 
for the implementation of this local modification computed from an electrostatic 
potential known on a spatio-temporal grid. The number of particles which escape 
the inner region defined by the barrier measures the efficiency of the control. 
We show that the control is robust by showing a significant reduction of radial 
transport, when applied with a reduced number of probes aligned on a circle. 



1. Introduction 

Control of chaotic transport is an important topic in many areas of physics with 
considerable research and technical applications. Its aim is to reduce chaos when it 
is harmful and enhance it when it is beneficial. A lot of efforts have been devoted to 
reduce chaotic and turbulent transport in laboratory plasma physics. Two examples 
are afforded by the control of stochasticity of magnetic field lines and the transport 
generated by plasma instabilities. In magnetically confined fusion plasmas, like 
tokamaks, the electrostatic turbulence associated with the E x B drift motion, is 
one of the main sources for the loss of energy confinement [2 [21 [3] and consequently 
constitutes a relevant obstacle to the attainment of plasma ignition. The development 
of control strategies able to induce a relevant change in transport properties by means 
of a small and localized perturbation constitutes a way to understand the complex 
behavior of laboratory plasmas. 

A convenient way to control these systems is by using electromagnetic fields. Since 
the plasma dynamics is very sensitive to electrostatic fluctuations, a control method 
for the E X B drift motion of test particles has been proposed by a small and apt 
modification of the electrostatic potential [H [5] . This method has been tested on a 
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Travelling- Wave Tube (TWT) for the reduction of chaotic transport [6,. It is worth 
noticing that for these methods, the relevant modifications were applied in all the 
plasma region. To step forward in the applicability of this control scheme to more 
complex plasma devices, we present a method and its practical algorithm to construct 
a transport barrier for the radial dynamics where the modification of the potential is 
only applied locally in a peripheral region of the plasma. Our strategy does not depend 
on the electrostatic potential at hand whether it is known analytically, numerically 
or experimentally. The geometry we are considering is a cylindrical spatial geometry 
which is shared by many linear devices like for example VINETA [7] , MIRABELLE [8] 
or CSDX [9]. Moreover, as in these devices the plasma temperature is not too high, 
plasma density and electrostatic potential are routinely measured by inserting probe 
arrays with large radial and poloidal extents. In that way it is possible to have a very 
accurate description of the space time behavior of the plasma [TOl [TTl [7] . In addition 
it is worth noticing that it is also possible to act on the plasma using electrodes 
directly in contact with the plasma, producing for example a mode locking between 
the injected sinusoidal signal and the plasma dynamics [121 [13]. Moreover, in the SOL 
region of tokamaks where the field lines are connected to the wall, a parallel current 
created by probes acts on perpendicular dynamics [H]. 

Section 2 deals with the analytical solution of the potential to create a radial transport 
barrier for test particles. In Sec. 3, we describe the numerical algorithm to compute 
this modified potential from the electrostatic potential data. Finally, Section 4 is 
dedicated to the computation and the results of the control of test particle dynamics 
in a mock electrostatic potential. 

2. Method 

2.1. Guiding- center dynamics 

The equations of motion of charged particles (of charge e and mass m) in 
electromagnetic fields (E = — VV — dA/dt and B = V x A) can be obtained in 
a Hamiltonian framework [151 IIS] where the Hamiltonian is 

2m 

where q is the position of the particle and p is its canonically conjugate momentum. 
In strong magnetic fields, particles have a fast gyration motion around the guiding- 
center [T7]. The relation between the position of the particle and its guiding center 
is: 

q = X + p, 

where x is the position of the guiding-center and p is the Larmor vector. In what 
follows, we restrict our study to the case where the magnetic field is constant and 
uniform: B = Bz. In this case, in order to reduce the system, we apply the generalized 
canonical transformation (q, p) ^ (x, u) given byx = q — zx(p — eA)/eB, and 
u = (p — eA)/m. The Hamiltonian becomes [T8l [T9]: 

and the Poisson bracket has pairs of conjugate coordinates: (z, Uz), {x, y) and {ux^Uy). 

- — z • (— X —\ + — z (— X —\ 
e5 \ 9x 9x y \ du du J 



Transport barrier for the radial diffusion due to the E x B drift 



3 



1 dF^ dG ^ dF ^ dG 



It is important to notice that the two coordinates {x and y) of the guiding center in 
the plane transverse to the magnetic field are conjugate, as well as the two velocities 
{ux and Uy) of the particle in this transverse plane. Nothing is changed along the 
direction of the field lines. In order to reduce the dimensionality of our system, we 
approximate the potential by 



which is a standard hypothesis for plasma devices with a strong magnetic field. In 
this way, we notice that the equations of motion for the position of the guiding center 
are decoupled from the ones of the velocity of the particle. The equations of motion 
are given by 

I dV 

" 

1 dV 

that is the velocity of the guiding centers is equal to the E x B drift velocity. 

We consider a cylindrical geometry as encountered in linear plasma devices. In 
this geometry, the polar coordinates are better suited. We perform the (non-canonical) 
change of coordinates x = rcosO and y = rs'mO. The Poisson bracket of the reduced 
system in the poloidal plane is changed into 

1 fdFdG dFdG 

with a reduced Hamiltonian H = eV{r^9^t)^ where Vir^O^t) = V{x^y^t). In what 
follows, we remove the tildas for simplicity of the notations. Hence, the equations of 
motion for the guiding centers are: 

The singularity at r = can be removed by reformulating the equations of motion 
in the canonically conjugate pair of variables (^,^) with ^ = However the 

singularity is in general not a problem in the numerical integration of test particles 
since the measured fiuctuating potentials vanish in the core region around r = 0. The 
important region where confinement is crucial is the peripheral (edge) region. 

2.2. Control term for a transport barrier 

We are looking at modifications of the Hamiltonian, i.e. of the electrostatic potential 
V, in order to create a radial transport barrier. We consider the following expression 
for the controlled potential Vc'. 

v,{r,e,t) = v{f{r,e,t),e,t), (3) 

where f{r,6,t) is the control function that has to be determined. We define the map 
T as: 

T:(r,^,t)^(/(r,^,t),^,t), (4) 
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with which the controhed potential is written as Vc = V oT. We impose that this new 
potential creates a transport barrier for the guiding centers located at r = R{0^t). 
This is equivalent to imposing that r = R{9^t) is invariant by the dynamics. More 
specifically, we have 

R = dtR{e,t)^OdeR{e,t), (5) 

where R and are given by Eqs. ([T]) and ([2]) for the potential Vc given by Eq. (|3]). 
Equation ([5j) translates into 

1 dV d . . dR I dV ^ 

where the derivatives of V are taken at r = f{R{0^t)^O^t). This equation has two 
unknown functions / and R to be determined. We assume that the time-average of V 
is zero for all ^ at a given radius ro- We obtain a solution for / and R assuming that 

f{R,0,t) = ro, (6) 

and 



R(e,t) = ^r^,-^r^{ro,e,t), (7) 

where F is a linear operator acting on any function (with vanishing mean value) of 
time as the pseudo-inverse of the derivative with respect to time, e.g., acting on a 
function (j){t) as 

There is still a lot of freedom for the choice of the control term, even if all these 
choices lead to the same equation for the transport barrier R given by Eq. ([7]). We 
solve Eq. (|6]) using the following control terms: 

hir,e,t)=r + MO,t), (9) 

f2{r,e,t)=r + e-^^--'>^"/^^"M0,t), (10) 

where the V^i(6>, t) are arbitrary functions of 6 and t to be determined by Eq. (j6j). The 
first expression /i for the control term implies that the potential has to be modified 
on the entire phase space. However it is desirable to only apply the modifications on 
a particular region of the transverse plane. In order to localize radially the control 
required for creating the transport barrier, a first idea is to localize the potential 
modifications in a neighborhood of the barrier R(6^t). The fact that this localization 
is varying in time makes this strategy challenging for an experimental realization. In 
order to address this issue we adopt a slightly different strategy: the control term 
/2 leads to a modification of the potential only in the neighborhood of a circle of 
radius ro which is more adapted to an experimental application since the control is 
envisaged through a set of probes usually aligned along a circle, or a predefined fixed 
geometry (see below). A Gaussian function around the radius ro is used to describe 
this localization and cr is related to the width at half maximum. Using Eqs. (j3j), ([6]) 
and ([7j), the controlled potentials associated with the control terms (j9j) and (p!Q|) are 
respectively: 

Vc{r,e,t) = V{r^r^-R,e,t), (11) 
K(r, e,t) = v(r^ (ro - R)e^^R-^o?-{r-ro?)/{2a') ^ ^ (12) 

where R{0^t) is defined by Eq. ([7]). 
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3. Numerical algorithm 

We investigate numerically the effect of the controlled potential Vc = V o T on 
the dynamics of guiding centers. We have developed our algorithm such that it 
is applicable for any potential known on a spatio-temporal grid as obtained from 
numerical codes or as measured experimentally. 

3.1. Computation of the controlled electrostatic potential 

The measurement of the electrostatic potential V is performed on a grid of x Mq x 
Mt points called the measurement grid. This is the input of our numerical algorithm. 
This measurement can be given by experimental measurements or numerical data. 
However the control term / is computed on a refined spatial grid of Nr x Ng x Nt 
points where A^^ > and Nq > Mq in order to minimize the error in the numerical 
scheme. The potential given by Eq. ^ is used in order to create a local barrier of 
radial transport. The controlled potential Vc is computed on the measurement grid 
or on the refined one if the original grid is too coarse-grained (in order to avoid fake 
dissipation due to an imprecise way of computing Hamilton's equations). The output 
of this computation is the modifications of the electrostatic potential given on the 
measurement grid. 

The derivatives with respect to periodic variables (e.g., 0) as well as the F-operator 
given by Eq. ([8]) are computed using the Fast Fourier Transform and for non-periodic 
variables (e.g., r) with a fourth order finite-difference scheme. 

As we can see from Eqs. (pTj) and ([12]), the controlled potential Vc is simply obtained 
from the original one F by a shift in the radial position r (which depends on and 
t). This means that for determining Vc on a given spatial grid, one has to recover 
by interpolation, the values of V on a deformed grid as schematically represented in 
Fig. [TJ In this way, the method does not resort to an expansion of the potential and is 
more easily applicable to electric potentials obtained numerically or experimentally. 

3.2. Computation of test particles 

We compute and compare the dynamics of test particles given by the equations 
of motion ([T]) and (|2]) obtained with the uncontrolled and controlled electrostatic 
potentials, V and Vc = V o T respectively. In order to integrate these equations 
of motion, we use a fourth order Runge-Kutta integrator. Even though this is not 
a symplectic scheme, we monitor the total energy E + V{x(t)^y{t)^t) where the 
dynamics of the additional variable £ is defined hy £ = —dV/dt^ so as to ensure 
that the numerical error is negligible. This also constitutes a way to check globally 
our numerical scheme. The dynamical equation for £ is integrated with the same 
numerical scheme as Eqs. ([T]) and (|2]). 

With a set of initial conditions inside the barrier r = R{0^t)^ we compute the 
dynamics of test particles separately for the uncontrolled and for the controlled 
electrostatic potentials. If there is radial transport, the particles which are initially 
inside the barrier may leave the inner region. When the control creates a radial 
transport barrier, the test particles are expected to remain inside this barrier at all 
times. 




Figure 1. The controlled potential Vc on the initial grid is obtained by 
interpolating the uncontrolled potential V (from the red grid) on the deformed 
blue grid obtained from the map T given by Eq. Q. 



4. Numerical results 

4'1- Computation with uncontrolled potential 

In order to validate our method with the computation of test particles, we use a 
toy analytic potential with a decay of the spatial modes in consistent with 
experimental measurements [21] |l], which corresponds to the inertial range of the 
drift wave turbulence spectrum. We use only one frequency given that the frequency 
spectrum does not introduce changes in the algorithm as shown in Ref. ^22^. This 
potential has been studied in Ref. and is given by: 

^ 1 

y(r, ^,t)=a ^ sin (nr cos ^ + mr sin ^ + (/p^^ — t) , (13) 

n,m — —N ^ ' 

+ < N'^ 
n ^ 0,m ^ 

where a is the amplitude of the potential and (pnm are random phases in order to model 
a turbulent potential with N = 25. The amplitude a is a parameter to modulate the 
importance of chaotic transport. Figure [2] is a representation of this potential at time 
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t = 0. We add a filter in a small circular region centered in r = to model the fact 
that the electrostatic fluctuations are weak at the center of the device. 




Figure 2. Potential V given by Eq. ([13]) at t = for a = 0.4. The initial 
condition of test particles positions are chosen randomly inside the barrier R{0, 0) 
which will be obtained for the controlled potential Vc (see Fig. [6])- 

In order to obtain an accurate numerical simulation, we choose the number of 
points Nr = 256 (in the radial position r), Nq = 512 (in the poloidal angle 6) 
and Nt = 64 (in time) for the electrostatic potential grid. The measurement of the 
electrostatic potential is performed with {Mr^Mg) = (128,256) and Mf = Nf. The 
total energy conservation is accurate up to 10~^ up to a time of integration of 100 
periods. The initial positions of test particles are chosen near the center (randomly 
inside the position of the barrier R{6^ 0) obtained later with the controlled electrostatic 
potential Vc, i.e. inside the interval [0.45; 2.55] where the minimum radius of the barrier 
will be approximately 2.6). The Poincare section (stroboscopic plot with the period of 
the potential) of test particles obtained with the uncontrolled electrostatic potential is 
shown in Fig. [H The green curve defined by r = R{6, t) at t = is added to show that 
guiding-centers can diffuse (without control) in all the available plane perpendicular to 
the magnetic field lines. More precisely, in the computation presented in Fig. [3l only 
12% of the trajectories remain inside the barrier after 100 periods of the potential. 
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Figure 3. Poincare section of 100 test particles for the potential V given by 
Eq. (|13|) for a = 0.4, up to t = 100 x 27r. The green curve is the position of the 
barrier R(0, 0) which will be obtained with the controlled potential Vc (see Fig. (6]). 

4.2. Computation with controlled potential 
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Figure 4. Initial grid (red dots) where the uncontrolled potential V is known. 
Modified grid (blue lines) defined by Eq. (|4]) where V is interpolated to obtain Vc 
given by Eq. (|11|) on the initial grid (red dots). The green curve is the position 
of the barrier fiuctuating near the fixed point ro = tt at time t = 0. We use the 
potential V given by Eq. (|13|) with a = 0.4. 

The control term and electrostatic potential are computed on the measurement 
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Figure 5. Difference Vc — V between Vc given by Eq. (|12|) with cr = 0.6 and V 
given by Eq. (|13|) at t = for a = 0.4. This modification is centered around the 
black circle with radius ro = tt. 



grid (Mr^Mg). An additional step is to interpolate the controlled potential on the 
initial grid {Nr^Ng) using a spline interpolation. We apply our control algorithm 
in order to confine test particles inside the barrier. The controlled potential given by 
Eq. (p!2]) is considered. It is computed numerically performing an interpolation method 
of V on the deformed grid explicitly defined hy Vc = V o T following the procedure 
sketched in Fig. [H For example, in Fig. |4] we represent the deformation of the grid 
associated with the computation of the controlled potential given by Eq. (pTj) . The 
difference Vc — V where V is given by Eq. (p!3|) and Vc is given by Eq. ([12]) is shown in 
Fig. [51 In this case, the maximum amplitude of the modification Vc — V is less than 
50% of the maximum amplitude of the potential V for a = 0.4. 

The efficiency of the control algorithm is estimated from the dynamics of test 
particles obtained by integrating Eqs. ([!]) and ([2]) with the controlled potential Vc. 
The Poincare section of test particles obtained with the controlled potential Vc given 
by Eq. ([T2]) is shown in Fig. [6l 

The position of the control barrier R{0^ 2/c7r) = R{0^ 0) is represented by the solid 
red curve. We notice that most of the particles do not cross the transport barrier 
r = RiO^t). A quantitative estimate of the efficiency of the control scheme is obtained 
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Figure 6. Poincare section of 100 test particles for the potential Vc given by 
Eq. (|12|) with a = 0.6 and ro = tt where V is given by Eq. ()13|) for a = 0.4 up to 
t = 100 X 27T. The red curve is the position of the barrier R(0, 0). 

by comparing the number of particles trapped inside the barrier for the controhed and 
uncontrolled potentials. The time evolution of the percentage of particles Nconfined 
which stay inside the barrier is shown in Fig. [71 For a = 0.4, up to t = 100 x 27r, there is 
12% of particles that remain confined within the prescribed region without the control 
(see Fig. [3]) in comparison to 90% for the controlled potential (see Fig. [6]). We notice 
that the control does not lead to 100% of trapped particles (the dashed blue curve 
in Fig. [7j) as expected from Eq. ([12]) since the numerical algorithm (performed on a 
grid) induces an approximation of the exact control. This is also to be expected from 
the discrepancy between the measured potential and the actual one which prevents 
the computation of a very accurate controlled potential. However the fact that a 
significant number of test particles are trapped inside the transport barrier shows 
that the proposed control is robust in the numerical implementation. This robustness 
allows us to take into account experimental constraints. Here we have considered 
that the controlled potential is applied through a set of probes. Two constraints 
were considered: First, the geometry with which the probes are introduced should 
be as simple as possible. Second, the number of probes necessary for an effective 
implementation of the controlled potential should be relatively small and compatible 
with what is currently applied in linear plasma devices. 

The first constraint is addressed by aligning the probes on a fixed circle of radius 
r = To (see Fig. [8]); the second one by using Mprobes = 32 probes and measuring its 
effects on the number of trapped particles. We assume that each probe creates an 
electrostatic potential with a Gaussian decrease of its infiuence in both the radial and 
poloidal direction. The potential created by the probe i at the position (r, 6) = (ro, 6>i) 
is given by 

V,{r,e,t) = A,(t)e-(^-^°)'/2^%-(^-^^)'/2^% (14) 
where Ai{t) is the amplitude of the electric potential of the probe i at time t and 
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Figure 7. Time evolution of the percentage of confined particles which stay inside 
the barrier i.e. with (r,0) such that r < aR(6, t) where a = 1.2. The parameters 
of the potential are the same as in Fig. (6] This picture shows the results without 
(in solid black curve) or with the control (in dashed blue curve) computed from 
2000 particles. The control obtained with 32 probes on a circle of radius r = ro 
is also shown (in dash-dotted red curve without measurement constraints or in 
dotted green curve with a measurement on 5 x 128 points). 

<j = 0.6 (respectively ae = 27i / Mprohes) is the spatial decrease of the electric potential 
in the radial (respectively poloidal) direction. 

If the potential introduced at r = ro is such that the implemented potential at 
r = R{0^ t) is as close as possible to the potential Vc — V depicted in Fig. [5] (by adjusting 
the amplitudes Ai{t)) using the ordinary least square method, as shown in Fig. [9l 
then our numerical results show that already with 32 probes, the control strategy 
is effective in reducing radial transport as shown in Fig. [3 More quantitatively, we 
found that at t = 100 x 27r, the number of trapped particles which is about 12% 
without the control is increased up to 29% with 32 probes without measurement 
constraints and up to 21% with 32 probes with a measurement on 5 x 128 grid points. 
In addition to robustness, an important advantage of the proposed strategy concerns 
the relatively small modifications of the potential. The mean average in angle and 
time t (denoted (•)) of the potential modifications on the radial position r = ro, 
i.e. (|y(ro,6>,t) — Vc{ro^O,t)\), is about 30% of the mean average of the potential 
(|y(ro, 0,t)\) for a = 0.4. This has to be compared with a simpler version of the control 
which consists in creating a potential barrier around r = ro such that Vc{r = r^^O^t) 
is small. In order to obtain the same quantity of trapped particles {Nconfined ^ 90% 
up to t = 100 X 27r), it would be necessary to cancel 90% of the electrostatic potential 
which has to be compared with 30% of the modifications introduced by the proposed 
strategy with the same parameters as in Fig. [6l 
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Figure 8. The perturbations of the electric potential created by 32 probes at 
t = 0. The parameters of the potential are the same as in Fig. (6] The gray dots 
show the position of the probes and the solid red curve is the position of the 
barrier R(0, 0). The 5 black crosses correspond to the position of one line of the 
measurement points which give the dotted green curve in Fig. [Jl 
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